%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Notes: This simulation combines linear supply curves for ethanol from
% grain and cellulosic biomass with perfectly elasticity gasoline supply and
% constant elasticity fuel demand and simulates emissions and carbon prices
% for the 2022 RFS and a LCFS and a cap and trade policy producing equivalent 
% emissions reductions
%
% Notation: RI - reference values, S - supply, D - demand, E -
% elasticities
%
% Version 35: Updated directories, added toggles for scenarios
% Version 34: Added tax into supply curve intercept
% Version 33: Added linear interpolation for the supply curves, JH added
% corn supply length back in
% Version 32: Put fast code in RFS loops (based on V30)
% Version 31: Put fuel matrix search loops back in, faster run (based on V29)
% Version 30: SPH removed fuel matrix search loops
% Version 29: Reporting PS including subsidies in EtOH stubsidy policy
% Version 28: CIe_corn = 0.80 adopted as preferred number
% Version 27: Updated LCFS and CT searches to two-stage search
% Version 26: Updated RFS calc's to deal with convergence issues (pick
% lambda upper after convergence)
% Version 25: Removed High/Low supply calculations
% Version 24: Major code clean up to remove vestigal functionality
% Version 23: Adds calculations for subsidies
% Version 22: Minor changes
% Version 21: Implemented correct RFS incentives in nested loops of 2*3=6
% Version 20: Simplified LCFS simulation to loop over lambda only, LCFS is
% implied and calculated from shadow price.  Added in VEETC and cellulosic
% ethanol subsidies, operated by a switch.  Changed oil price scenarios to
% gasoline price scenarios at retail prices of $2.25, $2.75, $3.25.
% Version 19: Reduced lambda lower starting value in LCFS simulation, added
% output statement for baseline values
% Version 18: Fixed LCFS simulation to loop over lambda (inner) and LCFS (outer)
% to hit RFS emissions
% Version 17 Notes: Same as version 14 but with msw split to msw and mse from 
% food waste.  MSW-food count as "advanced" not cellulosic under RFS.  Code
% has been streamlined for single model runs, i.e. no loops to gen. curves.
% Code has also been set up to sequentially solve RFS, then CT then LCFS
% that yield the same emissions.  Abatement and average abatement cost
% calculations have been added.  Functionality same as version 16 but code
% has been cleaned up a bit to remove vestigal lines.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear;   % Clear memory

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Simulation control parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

count_3_selected = 0;

% Set number of iterations
Max_iter = 20;  % Must be an even number
RFS_max_iter = 20;

% RFS Switch
% RFStoggle = 0 (skip RFS loops), RFStoggle = 1 (do RFS loops)
RFStoggle = 1;

% Innovation switch  (Note: this switch also shuts off the RFS loops)
% Innov_toggle = 0 (All technologies)  Innov_toggle = 1 (No cellulosic technologies)
Innov_toggle = 0;

% User Switch
user= 2;  % user=1 is Jon;  user=2 is Stephen


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Fuel Properties
BTUPGg = 116093;          % BTU/gal. LHV used to convert emissions/unit energy to gge
BTUPGe = 76330;         % BTU/gal. LHV don't curretly use this as ethanol supply curves are already in gge

% Elasticities
Egs = 2;                    % Price elasticity of gasoline supply (long-run)
Egged = -.5;                % Price elasticity of fuel btu demand (long-run)

% Transportation Fuel Market Initial Equilibrium (calibrated to 2008)
% Gasoline is United States total market 
PgRI = 2.75;    % Retail: (low = $2.25, base = $2.75, high = $3.25)
GGERI = 1.40E+11;  % Baseline = 1.40E+11 AEO 2022 value w/ new mpg "Old" = 1.34E+11 Baseline/Reference US gasoline consumption avg. (2000-2008) Product Supplied "Finished Motor Gasoline"


% Carbon intensities (ratios)
CIe_corn = 0.80; % Preferred = 0.80; High iLUC = 1.0; Waste Pass = 0.80; Old Corn = 0.90;
CIe_agres = 0.20; % Preferred = 0.20; High iLUC = 0.20; Waste Pass = 0.00; Old Corn = 0.20;
CIe_orch_vin = 0.20; % Preferred = 0.20; High iLUC = 0.20; Waste Pass = 0.00; Old Corn = 0.20;
CIe_forest = 0.20; % Preferred = 0.20; High iLUC = 0.20; Waste Pass = 0.00; Old Corn = 0.20;
CIe_hec = 0.25; % Preferred = 0.25; High iLUC = 0.40; Waste Pass = 0.25; Old Corn = 0.25
CIe_msw = 0.20; % Preferred = 0.20; High iLUC = 0.20; Waste Pass = 0.00; Old Corn = 0.20;
CIe_msw_food = 0.20; % Preferred = 0.20; High iLUC = 0.20; Waste Pass = 0.00; Old Corn = 0.20; 
CIg = 1.00; % Gasoline is the baseline fuel by definition for LCFS
  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Scenarios
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% scenario = 0 - Base case
% scenario = 1 - Low BAU fuel price
% scenario = 2 - High BAU fuel price
% scenario = 3 - High iLUC emissions scenario
% scenario = 4 - Waste Pass emissions scenario
% scenario = 5 - Old Corn emissions scenario
% scenario = 6 - Less elastic demand scenario
% scenario = 7 - More elastic demand scenario
% scenario = 8 - Corn price elasticity scenario

% Run one scenario
scenario = 0

if scenario == 1
PgRI = 2.25
end 
if scenario == 2
PgRI = 3.25
end 
if scenario == 3
CIe_corn = 1.0
CIe_hec = 0.40
end 
if scenario == 4
CIe_agres = 0.0
CIe_orch_vin = 0.0
CIe_forest = 0.0
CIe_msw = 0.0
CIe_msw_food = 0.0
end
if scenario == 5
CIe_corn = 0.9
end
if scenario == 6
Egged = -.3
end
if scenario == 7
Egged = -.7
end
if scenario == 8
corn_supply = csvread('/Users/Jon/Dropbox/LandUse/1 - Supply Curves/Individual Fuel csv/cornAdj.csv',1,0);
end


% Carbon Factors
CFg = 76.839+25.367;    % Gasoline carbon factor MMTCO2e/Quad (1 Quad = 10^15 Btu) CA GREET 1.8 Fuel plus upstream emissions (WTT sheet)
CFg = CFg*10^6*1000*1000/(10^15)*BTUPGg;    % Gasoline carbon factor in grams/gal (gge) 
CFe_corn = CIe_corn * CFg;      
CFe_agres = CIe_agres * CFg; 
CFe_orch_vin = CIe_orch_vin * CFg;
CFe_forest = CIe_forest * CFg;
CFe_hec = CIe_hec * CFg;
CFe_msw = CIe_msw * CFg;
CFe_msw_food = CIe_msw_food * CFg;

% Initialize fuel price to baseline fuel price (unregulated)
PggeRI = PgRI;   % Initialize gge price to initial gasoline price
Pgge = PgRI;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Supply Curve Matrices
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
if user==1
    if scenario ~= 8
    corn_supply = csvread('/Users/Jon/Dropbox/LandUse/1 - Supply Curves/Individual Fuel csv/corn.csv',1,0);
    end
    agres_supply = csvread('/Users/Jon/Dropbox/LandUse/1 - Supply Curves/Individual Fuel csv/agres.csv',1,0);
    orch_vin_supply = csvread('/Users/Jon/Dropbox/LandUse/1 - Supply Curves/Individual Fuel csv/orch_vin.csv',1,0);
    forest_supply = csvread('/Users/Jon/Dropbox/LandUse/1 - Supply Curves/Individual Fuel csv/forest.csv',1,0);
    msw_supply = csvread('/Users/Jon/Dropbox/LandUse/1 - Supply Curves/Individual Fuel csv/msw.csv',1,0);
    msw_food_supply = csvread('/Users/Jon/Dropbox/LandUse/1 - Supply Curves/Individual Fuel csv/msw_food.csv',1,0);
    hec_supply = csvread('/Users/Jon/Dropbox/LandUse/1 - Supply Curves/Individual Fuel csv/hec.csv',1,0);
end

if user==2
    if scenario ~= 8
    corn_supply = csvread('C:\AllStephen\My Dropbox\Jon n Chris\1 - Supply Curves\Individual Fuel csv\corn.csv',1,0);
    end
    agres_supply = csvread('C:\AllStephen\My Dropbox\Jon n Chris\1 - Supply Curves\Individual Fuel csv\agres.csv',1,0);
    orch_vin_supply = csvread('C:\AllStephen\My Dropbox\Jon n Chris\1 - Supply Curves\Individual Fuel csv\orch_vin.csv',1,0);
    forest_supply = csvread('C:\AllStephen\My Dropbox\Jon n Chris\1 - Supply Curves\Individual Fuel csv\forest.csv',1,0);
    msw_supply = csvread('C:\AllStephen\My Dropbox\Jon n Chris\1 - Supply Curves\Individual Fuel csv\msw.csv',1,0);
    msw_food_supply = csvread('C:\AllStephen\My Dropbox\Jon n Chris\1 - Supply Curves\Individual Fuel csv\\msw_food.csv',1,0);
    hec_supply = csvread('C:\AllStephen\My Dropbox\Jon n Chris\1 - Supply Curves\Individual Fuel csv\\hec.csv',1,0);
end

if Innov_toggle == 1
  agres_supply(:,2)=0;
    orch_vin_supply(:,2)=0;
    forest_supply(:,2)=0;
    msw_supply(:,2)=0;
    msw_food_supply(:,2)=0;
    hec_supply(:,2)=0;
RFStoggle = 0;
end

supply_length=length(hec_supply);  % Warning: this controls looping over the supply curves
supply_length_corn = length(corn_supply);

retail_tax_adj = 0.63 ; % Adjustement for taxes and retail markup = mean over 2000 - 2009
adj_mat = zeros(length(hec_supply),2); % Pick the length of arbitrary supply vector 
adj_mat(:,1)  = retail_tax_adj;

adj_mat_corn_adj = zeros(length(corn_supply),2); % Pick the length of arbitrary supply vector 
adj_mat_corn_adj(:,1)  = retail_tax_adj;
corn_supply = corn_supply + adj_mat_corn_adj;

% corn_supply = corn_supply + adj_mat;
agres_supply = agres_supply + adj_mat;
orch_vin_supply = orch_vin_supply + adj_mat;
forest_supply = forest_supply + adj_mat;
msw_supply = msw_supply + adj_mat;
msw_food_supply = msw_food_supply + adj_mat;
hec_supply = hec_supply + adj_mat;


%%%%%% begin calculate cost vectors
cost_corn_vec=zeros(supply_length_corn,1);
cost_corn_vec(1)=0.5*(retail_tax_adj+corn_supply(1,1))*corn_supply(1,2);
for i=2:supply_length_corn
    cost_corn_vec(i)=cost_corn_vec(i-1)+0.5*(corn_supply(i,1)+corn_supply(i-1,1))*(corn_supply(i,2)-corn_supply(i-1,2));
end

cost_agres_vec=zeros(supply_length,1);
cost_agres_vec(1)=0.5*(retail_tax_adj+agres_supply(1,1))*agres_supply(1,2);
for i=2:supply_length
    cost_agres_vec(i)=cost_agres_vec(i-1)+0.5*(agres_supply(i,1)+agres_supply(i-1,1))*(agres_supply(i,2)-agres_supply(i-1,2));
end

cost_orch_vin_vec=zeros(supply_length,1);
cost_orch_vin_vec(1)=0.5*(retail_tax_adj+orch_vin_supply(1,1))*orch_vin_supply(1,2);
for i=2:supply_length
    cost_orch_vin_vec(i)=cost_orch_vin_vec(i-1)+0.5*(orch_vin_supply(i,1)+orch_vin_supply(i-1,1))*(orch_vin_supply(i,2)-orch_vin_supply(i-1,2));
end

cost_forest_vec=zeros(supply_length,1);
cost_forest_vec(1)=0.5*(retail_tax_adj+forest_supply(1,1))*forest_supply(1,2);
for i=2:supply_length
    cost_forest_vec(i)=cost_forest_vec(i-1)+0.5*(forest_supply(i,1)+forest_supply(i-1,1))*(forest_supply(i,2)-forest_supply(i-1,2));
end

cost_msw_vec=zeros(supply_length,1);
cost_msw_vec(1)=0.5*(retail_tax_adj+msw_supply(1,1))*msw_supply(1,2);
for i=2:supply_length
    cost_msw_vec(i)=cost_msw_vec(i-1)+0.5*(msw_supply(i,1)+msw_supply(i-1,1))*(msw_supply(i,2)-msw_supply(i-1,2));
end

cost_msw_food_vec=zeros(supply_length,1);
cost_msw_food_vec(1)=0.5*(retail_tax_adj+msw_food_supply(1,1))*msw_food_supply(1,2);
for i=2:supply_length
    cost_msw_food_vec(i)=cost_msw_food_vec(i-1)+0.5*(msw_food_supply(i,1)+msw_food_supply(i-1,1))*(msw_food_supply(i,2)-msw_food_supply(i-1,2));
end

cost_hec_vec=zeros(supply_length,1);
cost_hec_vec(1)=0.5*(retail_tax_adj+hec_supply(1,1))*hec_supply(1,2);
for i=2:supply_length
    cost_hec_vec(i)=cost_hec_vec(i-1)+0.5*(hec_supply(i,1)+hec_supply(i-1,1))*(hec_supply(i,2)-hec_supply(i-1,2));
end
%%%%%% end calculate cost vectors
                   

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Solve for Ethanol Production in Unregulated Case (BAU)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

disp 'Solve for Ethanol Production in Unregulated Case (BAU)'

%Pgge = 20;      use for tester

PPriceVec=ones(7,1)*Pgge;

CalculateSupply;

ERI_corn=Se_corn;
cost_corn_RI= cost_corn;
ERI_agres=Se_agres;
cost_agres_RI= cost_agres;
ERI_orch_vin=Se_orch_vin;
cost_orch_vin_RI= cost_orch_vin;
ERI_forest=Se_forest;
cost_forest_RI= cost_forest;
ERI_msw=Se_msw;
cost_msw_RI= cost_msw;
ERI_msw_food=Se_msw_food;
cost_msw_food_RI= cost_msw_food;
ERI_hec=Se_hec;
cost_hec_RI= cost_hec;



% Calculate total ethanol supply in unregulated equilibrium
ERI = ERI_corn + ERI_agres + ERI_orch_vin + ERI_forest + ERI_hec + ERI_msw + ERI_msw_food;
    
% Gasoline Reference
cost_gas_RI = PggeRI * (GGERI - ERI); 

% Calculate total costs in unregulated case (BAU)
Cost_RI = cost_corn_RI+cost_agres_RI+cost_orch_vin_RI+cost_forest_RI+cost_hec_RI+cost_msw_RI+cost_msw_food_RI+cost_gas_RI;

GRI = GGERI - ERI;      % Gasoline share of unregulated fuel supply

% Baseline (BAU) carbon
Baseline_Carbon = (GRI*CFg/10^6)+(ERI_corn*CFe_corn/10^6)+(ERI_agres*CFe_agres/10^6)+(ERI_orch_vin*CFe_orch_vin/10^6)+(ERI_forest*CFe_forest/10^6)+(ERI_hec*CFe_hec/10^6)+(ERI_msw_food*CFe_msw_food/10^6)+(ERI_msw*CFe_msw_food/10^6); % CO2 eq. with LCFS - Result here is MTCO2e 

Baseline_Carbon_Intensity = ((GRI*CIg)+(ERI_corn*CIe_corn)+(ERI_agres*CIe_agres)+(ERI_orch_vin*CIe_orch_vin)+(ERI_forest*CIe_forest)+(ERI_hec*CIe_hec)+(ERI_msw_food*CIe_msw_food)+(ERI_msw*CIe_msw_food))/GGERI; 


% Consumer Surplus
% Area under demand curves
% Set lower bound for integral so that integral is finite
int_limit = (1/100000)*GGERI;  % lower limit as a percentage of gge consumption in the unregulated market
Dem_area_RI = PggeRI*(1/GGERI)^(1/Egged)*(1/(1/Egged+1))*GGERI^(1/Egged+1) - PggeRI*(1/GGERI)^(1/Egged)*(1/(1/Egged+1))*int_limit^(1/Egged+1);
CS_RI = PggeRI*(1/GGERI)^(1/Egged)*(1/(1/Egged+1))*GGERI^(1/Egged+1) - PggeRI*(1/GGERI)^(1/Egged)*(1/(1/Egged+1))*int_limit^(1/Egged+1) - (PggeRI*(GGERI-int_limit));


% Calculations
SS_RI = Dem_area_RI - Cost_RI;

% Producer Surplus
PS_RI = PggeRI*GGERI - Cost_RI;

PS_corn_RI = (PggeRI*ERI_corn-cost_corn_RI);
PS_agres_RI = (PggeRI*ERI_agres-cost_agres_RI);
PS_orch_vin_RI = (PggeRI*ERI_orch_vin-cost_orch_vin_RI);
PS_forest_RI = (PggeRI*ERI_forest-cost_forest_RI);
PS_msw_RI = (PggeRI*ERI_msw-cost_msw_RI);
PS_msw_food_RI = (PggeRI*ERI_msw_food-cost_msw_food_RI);
PS_hec_RI = (PggeRI*ERI_hec-cost_hec_RI);
PS_gas_RI = (PggeRI*GRI-cost_gas_RI);

% Save prive vector & CSPSVec
PPriceVec_BAU = PPriceVec;
CSPSVec_BAU= [CS_RI; PS_corn_RI; PS_agres_RI; PS_orch_vin_RI; PS_forest_RI; PS_msw_RI;PS_msw_food_RI;PS_hec_RI; PS_gas_RI; 9999999];


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input parameters -remaining modeling parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Create and initialize internal modeling variables
AFCI = 0;               % Initialize AFCI (firm's wgt. avg. carbon intensity) to 0
GRI = GGERI - ERI;      % Gasoline share of unregulated fuel supply
Pgge = PggeRI;




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Renewable Fuel Standard Simulation (RFS)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tic

if RFStoggle == 1
disp 'Renewable Fuel Standard Simulation (RFS)'

% Initialize internal variables
% AFCI = 0;               % Initialize AFCI to 0
Sg = GRI;               % Initialize gasoline supply
Se_corn = 0;            
Se_agres = 0;
Se_orch_vin = 0;        % Initialize ethanol supplies
Se_forest = 0;
Se_hec = 0;
Se_msw = 0;
Se_msw_food = 0;
Sgge = GGERI;           % Initialize supply of BTU's
% Pgge = (Pgge_upper+Pgge_lower)/2;
count1 = 0;
count2 = 0;
count3 = 0;     % Several counter variables
i = 1;
j = 1;

% Toggle = 1 forces simulation to overcomply with RFS
RFS_over = 0;

% Set RFS Constraints 
RFS_total = 36*10^9*(BTUPGe/BTUPGg); 
RFS_advanced = 21*10^9*(BTUPGe/BTUPGg);
RFS_cellulosic = 16*10^9*(BTUPGe/BTUPGg);

% Set RFS policy ratios based on BAU fuel production
% These are based on BAU gasoline production
RFS_total_ratio = RFS_total/GRI;
RFS_advanced_ratio = RFS_advanced/GRI;
RFS_cellulosic_ratio = RFS_cellulosic/GRI;

% Initialize Shadow Prices
lambda_corn_upper = 1;
lambda_corn_lower = 0;
lambda_corn = .5;
lambda_cell_upper = 1;
lambda_cell_lower = 0;
lambda_cell = .5;
lambda_adv_upper = 1;
lambda_adv_lower = 0;
lambda_adv = .5;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Need the following if ratios are ``updated'' by policy makers
% Comment out the following 22 lines and loop ends at bottom if 
% ratios are no updated (i.e. RFS volumetric targets are not met
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

count4 = 0;     % For corn (total) ratio
count5 = 0;     % For advanced ratio
count6 = 0;     % For cellulosic ratio

% Initialize total ratio
RFS_total_ratio_upper = 1;
RFS_total_ratio_lower = 0;
RFS_total_ratio = 0.5;

% Search for ratios such that volumetric targets are met
while count4 <= RFS_max_iter;

% Initialize advanced ratio
RFS_advanced_ratio_upper = 1;
RFS_advanced_ratio_lower = 0;
RFS_advanced_ratio = 0.5;
    
    while count5 <= RFS_max_iter;
        
    % Initialize cellulosic ratio
    RFS_cellulosic_ratio_upper = 1;
    RFS_cellulosic_ratio_lower = 0;
    RFS_cellulosic_ratio = 0.5;
        
        while count6 <= RFS_max_iter;

% End loops for ratio updating

                % Initialize total shadow price
                lambda_corn_upper = 1;
                lambda_corn_lower = 0;
                lambda_corn = .5;
        
                % Search for lambda's such that ratio's for corn, advanced and cellulosic fuels are met
                while count1 <= RFS_max_iter;
                    
                % Set lambda to upper bound to ensure
                % constraint satisfied (after convergence)
                if count1 == RFS_max_iter && RFS_over == 1
                    lambda_corn = lambda_corn_upper;
                end
                                     
                % Initialize advanced shadow price
                lambda_adv_upper = 1;
                lambda_adv_lower = 0;
                lambda_adv = .5;

                    while count2 <= RFS_max_iter;
                        
                    % Set lambda to upper bound to ensure
                    % constraint satisfied (after convergence)
                    if count2 == RFS_max_iter && RFS_over == 1
                        lambda_adv = lambda_adv_upper;
                    end
                      
                    % Initialize cellulosic shadow price
                    lambda_cell_upper = 1;
                    lambda_cell_lower = 0;
                    lambda_cell = .5;

                        while count3 <= RFS_max_iter;
                            
                            % Set lambda to upper bound to ensure
                            % constraint satisfied (after convergence)
                            if count3 == RFS_max_iter && RFS_over == 1
                                lambda_cell = lambda_cell_upper;
                            end

                            % Set fuel price based on RFS and MC gasoline

Pgge = PggeRI + lambda_cell*RFS_cellulosic_ratio + lambda_adv*RFS_advanced_ratio + lambda_corn*RFS_total_ratio;
PPriceVec=PPriceVec-PPriceVec+(Pgge+lambda_corn + lambda_cell + lambda_adv);
PPriceVec(1)=Pgge+lambda_corn;
PPriceVec(6)=Pgge+lambda_corn + lambda_adv;

%CalculateSupply;

%%%%%% begin calculating supply and costs for given PPriceVec
% Corn  PPriceVec(1)
count_3 = 1;
if (PPriceVec(1)<= 0)
    Se_corn = 0;
    count_3 = 9999; % Break out of loop
elseif (PPriceVec(1)< corn_supply(1,1))
    Se_corn = (PPriceVec(1)/ corn_supply(1,1)) * corn_supply(1,2);
    count_3 = 9999; % Break out of loop
end
if (PPriceVec(1) >= corn_supply(supply_length_corn,1))
    Se_corn = corn_supply(supply_length_corn,2);
    count_3 = 9999; % Break out of loop
end
while (count_3 <= supply_length_corn-1)
    if corn_supply(count_3+1,1) > PPriceVec(1)
        Se_corn = corn_supply(count_3,2)+ (PPriceVec(1)-corn_supply(count_3,1))/(corn_supply(count_3+1,1)-corn_supply(count_3,1))*(corn_supply(count_3+1,2)-corn_supply(count_3,2));
        count_3 = 9999; % Break out of loop
    end
    count_3 = count_3+1;
end

% Cellulosic Ethanol - Ag Residues PPriceVec(2)
count_3 = 1;
if (PPriceVec(2)<= 0)
    Se_agres = 0;
    count_3 = 9999; % Break out of loop
elseif (PPriceVec(2)< agres_supply(1,1))
    Se_agres = (PPriceVec(2)/ agres_supply(1,1)) * agres_supply(1,2);
    count_3 = 9999; % Break out of loop
end
if (PPriceVec(2) >= agres_supply(supply_length,1))
    Se_agres = agres_supply(supply_length,2);
    count_3 = 9999; % Break out of loop
end
while (count_3 <= supply_length-1)
    if agres_supply(count_3+1,1) > PPriceVec(2)
        Se_agres = agres_supply(count_3,2)+ (PPriceVec(2)-agres_supply(count_3,1))/(agres_supply(count_3+1,1)-agres_supply(count_3,1))*(agres_supply(count_3+1,2)-agres_supply(count_3,2));
        count_3 = 9999; % Break out of loop
    end
    count_3 = count_3+1;
end

% Cellulosic Ethanol - orchards and vineyard prunings  PPriceVec(3)
count_3 = 1;
if (PPriceVec(3)<= 0)
    Se_orch_vin = 0;
    count_3 = 9999; % Break out of loop
elseif (PPriceVec(3)< orch_vin_supply(1,1))
    Se_orch_vin = (PPriceVec(3)/ orch_vin_supply(1,1)) * orch_vin_supply(1,2);
    count_3 = 9999; % Break out of loop
end
if (PPriceVec(3) >= orch_vin_supply(supply_length,1))
    Se_orch_vin = orch_vin_supply(supply_length,2);
    count_3 = 9999; % Break out of loop
end
while (count_3 <= supply_length-1)
    if orch_vin_supply(count_3+1,1) > PPriceVec(3)
        Se_orch_vin = orch_vin_supply(count_3,2)+ (PPriceVec(3)-orch_vin_supply(count_3,1))/(orch_vin_supply(count_3+1,1)-orch_vin_supply(count_3,1))*(orch_vin_supply(count_3+1,2)-orch_vin_supply(count_3,2));
        count_3 = 9999; % Break out of loop
    end
    count_3 = count_3+1;
end

% Cellulosic Ethanol - forest thinning  PPriceVec(4)
count_3 = 1;
if (PPriceVec(4)<= 0)
    Se_forest = 0;
    count_3 = 9999; % Break out of loop
elseif (PPriceVec(4)< forest_supply(1,1))
    Se_forest = (PPriceVec(4)/ forest_supply(1,1)) * forest_supply(1,2);
    count_3 = 9999; % Break out of loop
end
if (PPriceVec(4) >= forest_supply(supply_length,1))
    Se_forest = forest_supply(supply_length,2);
    count_3 = 9999; % Break out of loop
end
while (count_3 <= supply_length-1)
    if forest_supply(count_3+1,1) > PPriceVec(4)
        Se_forest = forest_supply(count_3,2)+ (PPriceVec(4)-forest_supply(count_3,1))/(forest_supply(count_3+1,1)-forest_supply(count_3,1))*(forest_supply(count_3+1,2)-forest_supply(count_3,2));
        count_3 = 9999; % Break out of loop
    end
    count_3 = count_3+1;
end

% Cellulosic Ethanol - municipal solid waste PPriceVec(5)
count_3 = 1;
if (PPriceVec(5)<= 0)
    Se_msw = 0;
    count_3 = 9999; % Break out of loop
elseif (PPriceVec(5)< msw_supply(1,1))
    Se_msw = (PPriceVec(5)/ msw_supply(1,1)) * msw_supply(1,2);
    count_3 = 9999; % Break out of loop
end
if (PPriceVec(5) >= msw_supply(supply_length,1))
    Se_msw = msw_supply(supply_length,2);
    count_3 = 9999; % Break out of loop
end
while (count_3 <= supply_length-1)
    if msw_supply(count_3+1,1) > PPriceVec(5)
        Se_msw = msw_supply(count_3,2)+ (PPriceVec(5)-msw_supply(count_3,1))/(msw_supply(count_3+1,1)-msw_supply(count_3,1))*(msw_supply(count_3+1,2)-msw_supply(count_3,2));
        count_3 = 9999; % Break out of loop
    end
    count_3 = count_3+1;
end

% Advanced Ethanol - municipal solid waste from food  PPriceVec(6)
count_3 = 1;
if (PPriceVec(6)<= 0)
    Se_msw_food = 0;
    count_3 = 9999; % Break out of loop
elseif (PPriceVec(6)< msw_food_supply(1,1))
    Se_msw_food = (PPriceVec(6)/ msw_food_supply(1,1)) * msw_food_supply(1,2);
    count_3 = 9999; % Break out of loop
end
if (PPriceVec(6) >= msw_food_supply(supply_length,1))
    Se_msw_food = msw_food_supply(supply_length,2);
    count_3 = 9999; % Break out of loop
end
while (count_3 <= supply_length-1)
    if msw_food_supply(count_3+1,1) > PPriceVec(6)
        Se_msw_food = msw_food_supply(count_3,2)+ (PPriceVec(6)-msw_food_supply(count_3,1))/(msw_food_supply(count_3+1,1)-msw_food_supply(count_3,1))*(msw_food_supply(count_3+1,2)-msw_food_supply(count_3,2));
        count_3 = 9999; % Break out of loop
    end
    count_3 = count_3+1;
end

% Cellulosic Ethanol - herbaceous energy crops  PPriceVec(7)
count_3 = 1;
if (PPriceVec(7)<= 0)
    Se_hec = 0;
    count_3 = 9999; % Break out of loop
elseif (PPriceVec(7)< hec_supply(1,1))
    Se_hec = (PPriceVec(7)/ hec_supply(1,1)) * hec_supply(1,2);
    count_3 = 9999; % Break out of loop
end
if (PPriceVec(7) >= hec_supply(supply_length,1))
    Se_hec = hec_supply(supply_length,2);
    count_3 = 9999; % Break out of loop
end
while (count_3 <= supply_length-1)
    if hec_supply(count_3+1,1) > PPriceVec(7)
        Se_hec = hec_supply(count_3,2)+ (PPriceVec(7)-hec_supply(count_3,1))/(hec_supply(count_3+1,1)-hec_supply(count_3,1))*(hec_supply(count_3+1,2)-hec_supply(count_3,2));
        count_3 = 9999; % Break out of loop
    end
    count_3 = count_3+1;
end
%%%%%% end calculating supply and costs for given PPriceVec

  Dgge = GGERI*((Pgge/PggeRI)^Egged);      % Calculate demand for fuel at market price
  Se = Se_corn + Se_agres + Se_orch_vin + Se_forest + Se_hec + Se_msw + Se_msw_food;       % Calculate supply of ethanol
  Sg = Dgge - Se;      % Calculate supply of gasoline given perfectly elastic supply
  if Sg < 0,      disp   'Error supply of gas is negative (CT)'; err2,    end;
  Sgge = Sg + Se;     % Calculate supply of fuel (GGE's)

                            % Check for compliance with cellulosic ratio
                            Se_cellulosic = Se_agres + Se_forest + Se_orch_vin + Se_hec + Se_msw;
                            cellulosic_ratio = (Se_agres + Se_forest + Se_orch_vin + Se_hec + Se_msw)/Sg;

                            if cellulosic_ratio < RFS_cellulosic_ratio
                                lambda_cell_lower = lambda_cell;
                            else
                                lambda_cell_upper = lambda_cell;
                            end
                            lambda_cell = (lambda_cell_upper + lambda_cell_lower)/2;

                            count3 = count3 + 1;
                        end
                        % Set correct lambda
                        if RFS_over == 1
                        lambda_cell = lambda_cell_upper;
                        end
                        
                        % Reset counter
                        count3 = 0;

                            % Check for compliance with advanced fuel ratio
                            Se_oth_adv = Se_msw_food;
                            Se_advanced = Se_cellulosic + Se_oth_adv;
                            advanced_ratio = (Se_cellulosic + Se_oth_adv)/Sg;
                            if advanced_ratio < RFS_advanced_ratio
                                lambda_adv_lower = lambda_adv; 
                            else
                                lambda_adv_upper = lambda_adv;
                            end
                            lambda_adv = (lambda_adv_upper + lambda_adv_lower)/2;

                            count2 = count2 + 1;
                    end
                    % Set correct lambda
                    if RFS_over == 1
                    lambda_adv = lambda_adv_upper;
                    end
                    
                    % Reset counter
                    count2 = 0;    

                            % Check for compliance with total renewable fuel requirement
                            total_ratio = (Se_cellulosic + Se_oth_adv + Se_corn)/Sg;

                            if total_ratio < RFS_total_ratio
                                lambda_corn_lower = lambda_corn; 
                            else
                                lambda_corn_upper = lambda_corn;
                            end
                            lambda_corn = (lambda_corn_upper + lambda_corn_lower)/2;

                            count1 = count1 + 1;
                end
                % Set correct lambda
                if RFS_over == 1
                lambda_corn = lambda_corn_upper;
                end
                
                % Reset counter
                count1 = 0;            

                
  % Comment to "End ratio with updating" to do ratios w/out updating              
                
        % Check cellulosic volumes
        if Se_cellulosic < RFS_cellulosic
            RFS_cellulosic_ratio_lower = RFS_cellulosic_ratio;
        else
            RFS_cellulosic_ratio_upper = RFS_cellulosic_ratio;
        end 
        RFS_cellulosic_ratio = (RFS_cellulosic_ratio_upper + RFS_cellulosic_ratio_lower)/2;
        
        count6 = count6 + 1;
        
        end
        
        % Reset counter
        count6 = 0;
        
    % Check advanced volumes
    if Se_advanced < RFS_advanced
        RFS_advanced_ratio_lower = RFS_advanced_ratio;
    else
        RFS_advanced_ratio_upper = RFS_advanced_ratio;
    end 
    RFS_advanced_ratio = (RFS_advanced_ratio_upper + RFS_advanced_ratio_lower)/2;

    count5 = count5 + 1;
    
    end
    
    % Reset counter
    count5 = 0;
    
    
% Check total volumes
if Se < RFS_total
    RFS_total_ratio_lower = RFS_total_ratio;
else
    RFS_total_ratio_upper = RFS_total_ratio;
end 
RFS_total_ratio = (RFS_total_ratio_upper + RFS_total_ratio_lower)/2;

count4 = count4 + 1;

% Output simulation status
outer_loop = count4

RFS_loops(outer_loop,1) = outer_loop;
RFS_loops(outer_loop,2) = Se;
RFS_loops(outer_loop,3) = Se_corn;
RFS_loops(outer_loop,4) = Se_cellulosic;
RFS_loops(outer_loop,5) = RFS_total_ratio;
RFS_loops(outer_loop,6) = RFS_advanced_ratio;
RFS_loops(outer_loop,7) = RFS_cellulosic_ratio;
RFS_loops(outer_loop,8) = (Se_corn+Se_agres+Se_orch_vin+Se_forest+Se_msw+Se_msw_food+Se_hec)/Sg;
RFS_loops(outer_loop,9) = (Se_agres+Se_orch_vin+Se_forest+Se_msw+Se_msw_food+Se_hec)/Sg;
RFS_loops(outer_loop,10) = (Se_agres+Se_orch_vin+Se_forest+Se_msw+Se_hec)/Sg;
RFS_loops(outer_loop,11) = lambda_corn;
RFS_loops(outer_loop,12) = lambda_adv;
RFS_loops(outer_loop,13) = lambda_cell;

end
    
% Reset counter
count4 = 0;

% End ratio with updating

CalculateSupply;

% Carbon
C_RI = (GRI*CFg/10^6)+(ERI_corn*CFe_corn/10^6)+(ERI_agres*CFe_agres/10^6)+(ERI_orch_vin*CFe_orch_vin/10^6)+(ERI_forest*CFe_forest/10^6)+(ERI_hec*CFe_hec/10^6)+(ERI_msw_food*CFe_msw_food/10^6)+(ERI_msw*CFe_msw_food/10^6); % CO2 eq. with LCFS - Result here is MTCO2e 
C = (Sg*CFg/10^6)+(Se_corn*CFe_corn/10^6)+(Se_agres*CFe_agres/10^6)+(Se_orch_vin*CFe_orch_vin/10^6)+(Se_forest*CFe_forest/10^6)+(Se_hec*CFe_hec/10^6)+(Se_msw*CFe_msw/10^6)+(Se_msw_food*CFe_msw_food/10^6); % CO2 eq. with LCFS - Result here is MTCO2e

% Calculate change in private surplus and average abatement cost
% Total cost of production        

    % Gasoline
    cost_gas = PggeRI * Sg;

    % Total cost under policy
    Cost = cost_corn+cost_agres+cost_orch_vin+cost_forest+cost_hec+cost_msw+cost_msw_food+cost_gas;

    
% Area under demand curves
% Set lower bound for integral so that integral is finite
Dem_area = PggeRI*(1/GGERI)^(1/Egged)*(1/(1/Egged+1))*Dgge^(1/Egged+1) - PggeRI*(1/GGERI)^(1/Egged)*(1/(1/Egged+1))*int_limit^(1/Egged+1);


% Calculations
SS_RI = Dem_area_RI - Cost_RI;
SS = Dem_area - Cost;
chng_SS = SS-SS_RI;
chng_carbon = C-C_RI;
Avg_abatement_cost = chng_SS/chng_carbon;

% Consumer Surplus
CS_RI = PggeRI*(1/GGERI)^(1/Egged)*(1/(1/Egged+1))*GGERI^(1/Egged+1) - PggeRI*(1/GGERI)^(1/Egged)*(1/(1/Egged+1))*int_limit^(1/Egged+1) - (PggeRI*(GGERI-int_limit));
CS =    PggeRI*(1/GGERI)^(1/Egged)*(1/(1/Egged+1))* Dgge^(1/Egged+1) - PggeRI*(1/GGERI)^(1/Egged)*(1/(1/Egged+1))*int_limit^(1/Egged+1) - (Pgge*(Dgge-int_limit));
chng_CS = CS - CS_RI;

% Producer Surplus
PS_RI = PggeRI*GGERI - Cost_RI;
PS = Pgge*Sgge - Cost;
chng_PS = PS - PS_RI;

% Recalculate the ratios
 RFS_total_ratio=(Se_corn+Se_agres+Se_orch_vin+Se_forest+Se_msw+Se_msw_food+Se_hec)/Sg
 RFS_advanced_ratio=(Se_agres+Se_orch_vin+Se_forest+Se_msw+Se_msw_food+Se_hec)/Sg
 RFS_cellulosic_ratio=(Se_agres+Se_orch_vin+Se_forest+Se_msw+Se_hec)/Sg

% Producer Surplus by fuel type
chng_PS_corn_RFS = (Pgge*Se_corn-cost_corn)-(PggeRI*ERI_corn-cost_corn_RI)+(lambda_corn)*Se_corn
chng_PS_agres_RFS = (Pgge*Se_agres-cost_agres)-(PggeRI*ERI_agres-cost_agres_RI)+(lambda_corn+lambda_adv+lambda_cell)*Se_agres
chng_PS_orch_vin_RFS = (Pgge*Se_orch_vin-cost_orch_vin)-(PggeRI*ERI_orch_vin-cost_orch_vin_RI)+(lambda_corn+lambda_adv+lambda_cell)*Se_orch_vin
chng_PS_forest_RFS = (Pgge*Se_forest-cost_forest)-(PggeRI*ERI_forest-cost_forest_RI)+(lambda_corn+lambda_adv+lambda_cell)*Se_forest
chng_PS_msw_RFS = (Pgge*Se_msw-cost_msw)-(PggeRI*ERI_msw-cost_msw_RI)+(lambda_corn+lambda_adv+lambda_cell)*Se_msw
chng_PS_msw_food_RFS = (Pgge*Se_msw_food-cost_msw_food)-(PggeRI*ERI_msw_food-cost_msw_food_RI)+(lambda_corn+lambda_adv)*Se_msw_food
chng_PS_hec_RFS = (Pgge*Se_hec-cost_hec)-(PggeRI*ERI_hec-cost_hec_RI)+(lambda_corn+lambda_adv+lambda_cell)*Se_hec
chng_PS_gas_RFS = (Pgge*Sg-cost_gas)-(PggeRI*GRI-cost_gas_RI)-(lambda_cell*RFS_cellulosic_ratio + lambda_adv*RFS_advanced_ratio + lambda_corn*RFS_total_ratio)*Sg

chng_PS
chng_PS2 = chng_PS_corn_RFS + chng_PS_agres_RFS + chng_PS_orch_vin_RFS + chng_PS_forest_RFS + chng_PS_msw_RFS + chng_PS_msw_food_RFS + chng_PS_hec_RFS+ chng_PS_gas_RFS


% PS
% PS2=  (Pgge*Se_corn-cost_corn)+(lambda_corn)*Se_corn+ (Pgge*Se_agres-cost_agres)+(lambda_corn+lambda_adv+lambda_cell)*Se_agres+(Pgge*Se_orch_vin-cost_orch_vin)+(lambda_corn+lambda_adv+lambda_cell)*Se_orch_vin+(Pgge*Se_forest-cost_forest)+(lambda_corn+lambda_adv+lambda_cell)*Se_forest+(Pgge*Se_msw-cost_msw)+(lambda_corn+lambda_adv+lambda_cell)*Se_msw+(Pgge*Se_msw_food-cost_msw_food)+(lambda_corn+lambda_adv)*Se_msw_food+(Pgge*Se_hec-cost_hec)+(lambda_corn+lambda_adv+lambda_cell)*Se_hec+(Pgge*Sg-cost_gas)-(lambda_cell*RFS_cellulosic_ratio + lambda_adv*RFS_advanced_ratio + lambda_corn*RFS_total_ratio)*Sg
% 


% Save Parameters of Interest
Supply_Gas_RFS = Sg/10^9;         % Convert to billion gge
Supply_EtOH_RFS = Se/10^9;        % Convert to billion gge
Supply_GGE_RFS = Sgge/10^9;          % Convert to billion gge
Price_GGE_RFS = Pgge;
Shadow_Price_RFS_Total = lambda_corn;
Shadow_Price_RFS_Advanced = lambda_adv;
Shadow_Price_RFS_Cellulosic = lambda_cell;
Supply_EtOH_corn_RFS = Se_corn/10^9;          % Convert to billion gge
Supply_EtOH_agres_RFS = Se_agres/10^9;            % Convert to billion gge
Supply_EtOH_orch_vin_RFS = Se_orch_vin/10^9;          % Convert to billion gge
Supply_EtOH_forest_RFS = Se_forest/10^9;             % Convert to billion gge
Supply_EtOH_hec_RFS = Se_hec/10^9;            % Convert to billion gge
Supply_EtOH_msw_RFS = Se_msw/10^9;            % Convert to billion gge
Supply_EtOH_msw_food_RFS = Se_msw_food/10^9;            % Convert to billion gge
Supply_EtOH_advanced = Se_advanced/10^9;            % Convert to billion gge
Supply_EtOH_cellulosic = Se_cellulosic/10^9;            % Convert to billion gge
Baseline_Carbon_RFS = C_RI/10^6;        % In million metric tons
Carbon_RFS= C/10^6;                     % In million metric tons
Change_CS_RFS = chng_CS/10^9;           % In billions of $
Change_PS_RFS = chng_PS/10^9;           % In billions of $
Abatement_Cost_RFS = chng_SS/10^9;           % In billions of $
Avg_Abatement_cost_RFS = Avg_abatement_cost;  % In $/MTCO2e

% Save prive vector & CSPSVec
PPriceVec_RFS = PPriceVec;
CSPSVec_RFS= [chng_CS;chng_PS_corn_RFS;chng_PS_agres_RFS;chng_PS_orch_vin_RFS;chng_PS_forest_RFS;chng_PS_msw_RFS;chng_PS_msw_food_RFS;chng_PS_hec_RFS;chng_PS_gas_RFS;9999999];

end       % end "if RFStoggle==1"

toc


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Carbon Trading (CAT = Historical Baseline LCFS Simulation)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp 'Carbon Trading (CAT = Historical Baseline LCFS Simulation)'


% Initialize internal variables
% AFCI = 0;               % Initialize AFCI to 0
Sg = GRI;               % Initialize gasoline supply
Se_corn = 0;            
Se_agres = 0;
Se_orch_vin = 0;        % Initialize ethanol supplies
Se_forest = 0;
Se_hec = 0;
Se_msw = 0;
Se_msw_food = 0;
Sgge = GGERI;           % Initialize supply of BTU's
% Pgge_upper = 25;        
% Pgge_lower = 0;         % Set range for gge price search
Pgge = PggeRI;
lambda_upper = 3;     % If change, change in "check lambda_upper" below
lambda_lower = 0;     % Set initial range lambda search 
lambda = (lambda_upper + lambda_lower)/2;

if RFStoggle==1; carbon_target = Carbon_RFS*10^6; end;  % Convert carbon from MMT CO2e to MT CO2e and set target
if RFStoggle==0; carbon_target = 1453.7*10^6; end;  % Convert carbon from MMT CO2e to MT CO2e and set target

if Baseline_Carbon< carbon_target, disp 'Error: CnT Carbon target is met for lambda=0'; err1, end;
%%  Note:   lambda_upper is checked  after looping.  See "check lambda-upper"

for count1=1:Max_iter;
    % Set fuel price to MC of gasoline plus carbon tax/charge
    Pgge = PggeRI + lambda*(CIg);
    PPriceVec=[Pgge-lambda*CIe_corn;Pgge-lambda*CIe_agres;Pgge-lambda*CIe_orch_vin;Pgge-lambda*CIe_forest;Pgge-lambda*CIe_msw;Pgge-lambda*CIe_msw_food;Pgge-lambda*CIe_hec];

CalculateSupply;

  %
  Dgge = GGERI*((Pgge/PggeRI)^Egged);      % Calculate demand for fuel at market price
  Se = Se_corn + Se_agres + Se_orch_vin + Se_forest + Se_hec + Se_msw + Se_msw_food;       % Calculate supply of ethanol
  Sg = Dgge - Se;      % Calculate supply of gasoline given perfectly elastic supply
   if Sg < 0,  disp     'Error supply of gas is negative (CT)'; err4,     end;
  Sgge = Sg + Se;     % Calculate supply of fuel (GGE's)
  Carbon_CT = (Sg*CFg/10^6)+(Se_corn*CFe_corn/10^6)+(Se_agres*CFe_agres/10^6)+(Se_orch_vin*CFe_orch_vin/10^6)+(Se_forest*CFe_forest/10^6)+(Se_hec*CFe_hec/10^6)+(Se_msw*CFe_msw/10^6)+(Se_msw_food*CFe_msw_food/10^6);

    % Set new bounds for lambda
    if Carbon_CT > carbon_target
        lambda_lower = lambda;
    else
        lambda_upper = lambda;
    end
    lambda = (lambda_upper+lambda_lower)/2;
end

% check lambda_upper
if lambda_upper > 3-.0001,  disp 'Error supply of gas is negative (CT)'; err4,end

% Final calculations
lambda = lambda_upper;
Pgge = PggeRI + lambda*(CIg);
PPriceVec=[Pgge-lambda*CIe_corn;Pgge-lambda*CIe_agres;Pgge-lambda*CIe_orch_vin;Pgge-lambda*CIe_forest;Pgge-lambda*CIe_msw;Pgge-lambda*CIe_msw_food;Pgge-lambda*CIe_hec];

CalculateSupply;

%
  Dgge = GGERI*((Pgge/PggeRI)^Egged);      % Calculate demand for fuel at market price
  Se = Se_corn + Se_agres + Se_orch_vin + Se_forest + Se_hec + Se_msw + Se_msw_food;       % Calculate supply of ethanol
  Sg = Dgge - Se;      % Calculate supply of gasoline given perfectly elastic supply
   if Sg < 0,  disp     'Error supply of gas is negative (CT)'; err4,     end;
  Sgge = Sg + Se;     % Calculate supply of fuel (GGE's)
  Carbon_CT = (Sg*CFg/10^6)+(Se_corn*CFe_corn/10^6)+(Se_agres*CFe_agres/10^6)+(Se_orch_vin*CFe_orch_vin/10^6)+(Se_forest*CFe_forest/10^6)+(Se_hec*CFe_hec/10^6)+(Se_msw*CFe_msw/10^6)+(Se_msw_food*CFe_msw_food/10^6);

        
% Carbon  (sph:  Are these needed???)
C_RI = (GRI*CFg/10^6)+(ERI_corn*CFe_corn/10^6)+(ERI_agres*CFe_agres/10^6)+(ERI_orch_vin*CFe_orch_vin/10^6)+(ERI_forest*CFe_forest/10^6)+(ERI_hec*CFe_hec/10^6)+(ERI_msw*CFe_msw/10^6)+(ERI_msw_food*CFe_msw_food/10^6); % CO2 eq. with LCFS - Result here is MMTCO@e 
C = (Sg*CFg/10^6)+(Se_corn*CFe_corn/10^6)+(Se_agres*CFe_agres/10^6)+(Se_orch_vin*CFe_orch_vin/10^6)+(Se_forest*CFe_forest/10^6)+(Se_hec*CFe_hec/10^6)+(Se_msw*CFe_msw/10^6)+(Se_msw_food*CFe_msw_food/10^6); % CO2 eq. with LCFS - Result here is MTCO@e


% Calculate change in private surplus and average abatement cost

    % Gasoline
    cost_gas = PggeRI * Sg;
    
    % Calculate total production costs under policy
    Cost = cost_corn+cost_agres+cost_orch_vin+cost_forest+cost_hec+cost_msw+cost_msw_food+cost_gas;

% Area under demand curves
% Set lower bound for integral (so that integral is finite
Dem_area    = PggeRI*(1/GGERI)^(1/Egged)*(1/(1/Egged+1))*Dgge^(1/Egged+1) - PggeRI*(1/GGERI)^(1/Egged)*(1/(1/Egged+1))*int_limit^(1/Egged+1);

% Calculations
SS_RI = Dem_area_RI - Cost_RI;
SS = Dem_area - Cost;
chng_SS = SS-SS_RI;
chng_carbon = C-C_RI;
Avg_abatement_cost = chng_SS/chng_carbon;

% Permit/Allowance Revenue
CT_Revenue = lambda*(CIg*Sg + CIe_corn*Se_corn + CIe_agres*Se_agres + CIe_orch_vin*Se_orch_vin + CIe_forest*Se_forest + CIe_hec*Se_hec + CIe_msw*Se_msw + CIe_msw_food*Se_msw_food);

% Consumer Surplus
CS_RI = PggeRI*(1/GGERI)^(1/Egged)*(1/(1/Egged+1))*GGERI^(1/Egged+1) - PggeRI*(1/GGERI)^(1/Egged)*(1/(1/Egged+1))*int_limit^(1/Egged+1) - (PggeRI*(GGERI-int_limit));
CS =  PggeRI*(1/GGERI)^(1/Egged)*(1/(1/Egged+1))*Dgge^(1/Egged+1) - PggeRI*(1/GGERI)^(1/Egged)*(1/(1/Egged+1))*int_limit^(1/Egged+1) - (Pgge*(Dgge-int_limit));
chng_CS = CS - CS_RI;

% Producer Surplus
PS_RI = PggeRI*GGERI - Cost_RI;
PS = Pgge*Sgge - Cost;
PS = PS - CT_Revenue;
chng_PS = PS - PS_RI;

% Producer Surplus by fuel type
chng_PS_corn_CT = (Pgge*Se_corn-cost_corn)-(PggeRI*ERI_corn-cost_corn_RI)-lambda*CIe_corn*Se_corn;
chng_PS_agres_CT = (Pgge*Se_agres-cost_agres)-(PggeRI*ERI_agres-cost_agres_RI)-lambda*CIe_agres*Se_agres;
chng_PS_orch_vin_CT = (Pgge*Se_orch_vin-cost_orch_vin)-(PggeRI*ERI_orch_vin-cost_orch_vin_RI)-lambda*CIe_orch_vin*Se_orch_vin;
chng_PS_forest_CT = (Pgge*Se_forest-cost_forest)-(PggeRI*ERI_forest-cost_forest_RI)-lambda*CIe_forest*Se_forest;
chng_PS_msw_CT = (Pgge*Se_msw-cost_msw)-(PggeRI*ERI_msw-cost_msw_RI)-lambda*CIe_msw*Se_msw;
chng_PS_msw_food_CT = (Pgge*Se_msw_food-cost_msw_food)-(PggeRI*ERI_msw_food-cost_msw_food_RI)-lambda*CIe_msw_food*Se_msw_food;
chng_PS_hec_CT = (Pgge*Se_hec-cost_hec)-(PggeRI*ERI_hec-cost_hec_RI)-lambda*CIe_hec*Se_hec;
chng_PS_gas_CT = (Pgge*Sg-cost_gas)-(PggeRI*GRI-cost_gas_RI)-lambda*CIg*Sg;
chng_PS2 = chng_PS_corn_CT + chng_PS_agres_CT + chng_PS_orch_vin_CT + chng_PS_forest_CT + chng_PS_msw_CT + chng_PS_msw_food_CT + chng_PS_hec_CT+ chng_PS_gas_CT;

% Save data for combined graphs
Supply_Gas_CT = Sg/10^9;         % Convert to billion gge
Supply_EtOH_CT = Se/10^9;        % Convert to billion gge
Supply_GGE_CT = Sgge/10^9;          % Convert to billion gge
Price_GGE_CT = Pgge;
Shadow_Price_CT = lambda/(CFg/(10^6));        % Units of $/MT
Supply_EtOH_corn_CT = Se_corn/10^9;          % Convert to billion gge
Supply_EtOH_agres_CT = Se_agres/10^9;            % Convert to billion gge
Supply_EtOH_orch_vin_CT = Se_orch_vin/10^9;          % Convert to billion gge
Supply_EtOH_forest_CT = Se_forest/10^9;             % Convert to billion gge
Supply_EtOH_hec_CT = Se_hec/10^9;            % Convert to billion gge
Supply_EtOH_msw_CT = Se_msw/10^9;            % Convert to billion gge
Supply_EtOH_msw_food_CT = Se_msw_food/10^9;            % Convert to billion gge
Baseline_Carbon_CT = C_RI/10^6; % Convert to million metric tons
Carbon_CT= C/10^6;                  % Convert to million metric tons
Change_CS_CT = chng_CS/10^9;           % In billions of $
Change_PS_CT = chng_PS/10^9;           % In billions of $
CT_Revenue_CT = CT_Revenue/10^9;
Abatement_Cost_CT = chng_SS/10^9;           % In billions of $
Avg_Abatement_cost_CT = Avg_abatement_cost;  % In $/MTCO2e

% Save prive vector & CSPSVec
PPriceVec_CT = PPriceVec;
CSPSVec_CT= [chng_CS;chng_PS_corn_CT;chng_PS_agres_CT;chng_PS_orch_vin_CT;chng_PS_forest_CT;chng_PS_msw_CT;chng_PS_msw_food_CT;chng_PS_hec_CT;chng_PS_gas_CT;CT_Revenue];



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Energy Based LCFS Simulation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp 'Energy Based LCFS Simulation'

Sg = GRI;               % Initialize gasoline supply
Se_corn = 0;            
Se_agres = 0;
Se_orch_vin = 0;        % Initialize ethanol supplies
Se_forest = 0;
Se_hec = 0;
Se_msw = 0;
Se_msw_food = 0;
Sgge = GGERI;           % Initialize supply of BTU's
LCFS_upper = Baseline_Carbon_Intensity;
LCFS_lower = 0.75;       % If changed, also change 0.75 in hard coding below  See "Check LCFS_lower and LCFS_upper"
if Innov_toggle == 0
    upper_limit = 50;   % Set initial range lambda search  (set at 50 for main runs, 250 for Innovation)
else
    upper_limit = 250;
end
lower_limit = 0;
count1 = 0;

% check LCFS_upper
if Baseline_Carbon< carbon_target, disp 'Error: LCFS Carbon target is met for LCFS_upper'; err5, end;
% LCFS_lower and LCFS_upper are checked below.  See "Check LCFS_lower and LCFS_upper"

% Main LCFS loop
LCFS = (LCFS_upper + LCFS_lower)/2;
for count1=1:Max_iter
    lambda_upper = upper_limit;
    lambda_lower = lower_limit;
    lambda = (lambda_upper+lambda_lower)/2;
    %
    for count2=1:Max_iter
        Pgge = PggeRI + lambda*(CIg-LCFS);
        PPriceVec=[Pgge-lambda*(CIe_corn-LCFS);Pgge-lambda*(CIe_agres-LCFS);Pgge-lambda*(CIe_orch_vin-LCFS);Pgge-lambda*(CIe_forest-LCFS);Pgge-lambda*(CIe_msw-LCFS);Pgge-lambda*(CIe_msw_food-LCFS);Pgge-lambda*(CIe_hec-LCFS)];

        CalculateSupply;
%
  Dgge = GGERI*((Pgge/PggeRI)^Egged);      % Calculate demand for fuel at market price
  Se = Se_corn + Se_agres + Se_orch_vin + Se_forest + Se_hec + Se_msw + Se_msw_food;       % Calculate supply of ethanol
  Sg = Dgge - Se;      % Calculate supply of gasoline given perfectly elastic supply
    if Sg < 0,  disp     'Error supply of gas is negative (CT)'; lambda; LCFS; err5,    end;
  Sgge = Sg + Se;     % Calculate supply of fuel (GGE's)
  test_CI = (CIg*Sg + CIe_corn*Se_corn + CIe_agres*Se_agres + CIe_orch_vin*Se_orch_vin + CIe_forest*Se_forest + CIe_hec*Se_hec + CIe_msw*Se_msw + CIe_msw_food*Se_msw_food)/(Sg + Se_corn + Se_agres + Se_orch_vin + Se_forest + Se_hec + Se_msw + Se_msw_food);
%
  if test_CI<LCFS, lambda_upper=lambda; else lambda_lower=lambda; end; 
  lambda = (lambda_upper+lambda_lower)/2;
  if lambda<lower_limit+.0001 || lambda>upper_limit-.0001, disp  'Error: LCFS not properly converged'; err, end;
end;
  Carbon_LCFS = (Sg*CFg/10^6)+(Se_corn*CFe_corn/10^6)+(Se_agres*CFe_agres/10^6)+(Se_orch_vin*CFe_orch_vin/10^6)+(Se_forest*CFe_forest/10^6)+(Se_hec*CFe_hec/10^6)+(Se_msw*CFe_msw/10^6)+(Se_msw_food*CFe_msw_food/10^6);
if Carbon_LCFS < carbon_target, LCFS_lower=LCFS; else LCFS_upper=LCFS; end;
LCFS = (LCFS_upper + LCFS_lower)/2;
end;

% Check LCFS_lower and LCFS_upper 
if LCFS<0.75+.0001 || LCFS>Baseline_Carbon_Intensity-.0001, disp  'Error: LCFS not properly converged hitting carbon target'; err,  end;


%%% Final calculations  
LCFS=LCFS_lower;   % We need to recalculate lambda here
lambda_upper = upper_limit;
lambda_lower = lower_limit;     
lambda = (lambda_upper+lambda_lower)/2;
for count1=1:Max_iter
  Pgge = PggeRI + lambda*(CIg-LCFS);
  PPriceVec=[Pgge-lambda*(CIe_corn-LCFS);Pgge-lambda*(CIe_agres-LCFS);Pgge-lambda*(CIe_orch_vin-LCFS);Pgge-lambda*(CIe_forest-LCFS);Pgge-lambda*(CIe_msw-LCFS);Pgge-lambda*(CIe_msw_food-LCFS);Pgge-lambda*(CIe_hec-LCFS)];
  CalculateSupply;
%
  Dgge = GGERI*((Pgge/PggeRI)^Egged);      % Calculate demand for fuel at market price
  Se = Se_corn + Se_agres + Se_orch_vin + Se_forest + Se_hec + Se_msw + Se_msw_food;       % Calculate supply of ethanol
  Sg = Dgge - Se;      % Calculate supply of gasoline given perfectly elastic supply
    if Sg < 0,  disp     'Error supply of gas is negative (CT)';  err5,    end;
  Sgge = Sg + Se;     % Calculate supply of fuel (GGE's)
  test_CI = (CIg*Sg + CIe_corn*Se_corn + CIe_agres*Se_agres + CIe_orch_vin*Se_orch_vin + CIe_forest*Se_forest + CIe_hec*Se_hec + CIe_msw*Se_msw + CIe_msw_food*Se_msw_food)/(Sg + Se_corn + Se_agres + Se_orch_vin + Se_forest + Se_hec + Se_msw + Se_msw_food);

  if test_CI>LCFS_lower, lambda_lower=lambda; else lambda_upper=lambda; end; 
  lambda = (lambda_upper+lambda_lower)/2;
end;

lambda=lambda_upper;
  Pgge = PggeRI + lambda*(CIg-LCFS);
  PPriceVec=[Pgge-lambda*(CIe_corn-LCFS);Pgge-lambda*(CIe_agres-LCFS);Pgge-lambda*(CIe_orch_vin-LCFS);Pgge-lambda*(CIe_forest-LCFS);Pgge-lambda*(CIe_msw-LCFS);Pgge-lambda*(CIe_msw_food-LCFS);Pgge-lambda*(CIe_hec-LCFS)];
  CalculateSupply;
%
  Dgge = GGERI*((Pgge/PggeRI)^Egged);      % Calculate demand for fuel at market price
  Se = Se_corn + Se_agres + Se_orch_vin + Se_forest + Se_hec + Se_msw + Se_msw_food;       % Calculate supply of ethanol
  Sg = Dgge - Se;      % Calculate supply of gasoline given perfectly elastic supply
    if Sg < 0,  disp     'Error supply of gas is negative (CT)'; lambda; LCFS; err5,    end;
  Sgge = Sg + Se;     % Calculate supply of fuel (GGE's)
  test_CI = (CIg*Sg + CIe_corn*Se_corn + CIe_agres*Se_agres + CIe_orch_vin*Se_orch_vin + CIe_forest*Se_forest + CIe_hec*Se_hec + CIe_msw*Se_msw + CIe_msw_food*Se_msw_food)/(Sg + Se_corn + Se_agres + Se_orch_vin + Se_forest + Se_hec + Se_msw + Se_msw_food);
%
  Carbon_LCFS = (Sg*CFg/10^6)+(Se_corn*CFe_corn/10^6)+(Se_agres*CFe_agres/10^6)+(Se_orch_vin*CFe_orch_vin/10^6)+(Se_forest*CFe_forest/10^6)+(Se_hec*CFe_hec/10^6)+(Se_msw*CFe_msw/10^6)+(Se_msw_food*CFe_msw_food/10^6);



% Carbon in Metric Tons CO2e
C_RI = (GRI*CFg/10^6)+(ERI_corn*CFe_corn/10^6)+(ERI_agres*CFe_agres/10^6)+(ERI_orch_vin*CFe_orch_vin/10^6)+(ERI_forest*CFe_forest/10^6)+(ERI_hec*CFe_hec/10^6)+(ERI_msw*CFe_msw/10^6)+(ERI_msw_food*CFe_msw_food/10^6); % CO2 eq. with LCFS - Result here is MMTCO@e 
C = (Sg*CFg/10^6)+(Se_corn*CFe_corn/10^6)+(Se_agres*CFe_agres/10^6)+(Se_orch_vin*CFe_orch_vin/10^6)+(Se_forest*CFe_forest/10^6)+(Se_hec*CFe_hec/10^6)+(Se_msw*CFe_msw/10^6)+(Se_msw_food*CFe_msw_food/10^6); % CO2 eq. with LCFS - Result here is MMTCO@e


% Gasoline
cost_gas = PggeRI * Sg;

% Calculate total production costs under policy
Cost = cost_corn+cost_agres+cost_orch_vin+cost_forest+cost_hec+cost_msw+cost_msw_food+cost_gas;


% Area under demand curves
% Set lower bound for integral (so that integral is finite
Dem_area = PggeRI*(1/GGERI)^(1/Egged)*(1/(1/Egged+1))*Dgge^(1/Egged+1) - PggeRI*(1/GGERI)^(1/Egged)*(1/(1/Egged+1))*int_limit^(1/Egged+1);

% Calculations
SS_RI = Dem_area_RI - Cost_RI;
SS = Dem_area - Cost;
chng_SS = SS-SS_RI;
chng_carbon = C-C_RI;
Avg_abatement_cost = chng_SS/chng_carbon;

% Consumer Surplus
CS_RI = PggeRI*(1/GGERI)^(1/Egged)*(1/(1/Egged+1))*GGERI^(1/Egged+1) - PggeRI*(1/GGERI)^(1/Egged)*(1/(1/Egged+1))*int_limit^(1/Egged+1) - (PggeRI*(GGERI-int_limit));
CS =  PggeRI*(1/GGERI)^(1/Egged)*(1/(1/Egged+1))*Dgge^(1/Egged+1) - PggeRI*(1/GGERI)^(1/Egged)*(1/(1/Egged+1))*int_limit^(1/Egged+1) - (Pgge*(Dgge-int_limit));
chng_CS = CS - CS_RI;

% Producer Surplus
PS_RI = PggeRI*GGERI - Cost_RI;
PS = Pgge*Sgge - Cost;
chng_PS = PS - PS_RI;

% Producer Surplus by fuel type
chng_PS_corn_LCFS  = (Pgge*Se_corn-cost_corn)-(PggeRI*ERI_corn-cost_corn_RI)-lambda*(CIe_corn-LCFS)*Se_corn;
chng_PS_agres_LCFS = (Pgge*Se_agres-cost_agres)-(PggeRI*ERI_agres-cost_agres_RI)-lambda*(CIe_agres-LCFS)*Se_agres;
chng_PS_orch_vin_LCFS = (Pgge*Se_orch_vin-cost_orch_vin)-(PggeRI*ERI_orch_vin-cost_orch_vin_RI)-lambda*(CIe_orch_vin-LCFS)*Se_orch_vin;
chng_PS_forest_LCFS = (Pgge*Se_forest-cost_forest)-(PggeRI*ERI_forest-cost_forest_RI)-lambda*(CIe_forest-LCFS)*Se_forest;
chng_PS_msw_LCFS = (Pgge*Se_msw-cost_msw)-(PggeRI*ERI_msw-cost_msw_RI)-lambda*(CIe_msw-LCFS)*Se_msw;
chng_PS_msw_food_LCFS = (Pgge*Se_msw_food-cost_msw_food)-(PggeRI*ERI_msw_food-cost_msw_food_RI)-lambda*(CIe_msw_food-LCFS)*Se_msw_food;
chng_PS_hec_LCFS = (Pgge*Se_hec-cost_hec)-(PggeRI*ERI_hec-cost_hec_RI)-lambda*(CIe_hec-LCFS)*Se_hec;
chng_PS_gas_LCFS = (Pgge*Sg-cost_gas)-(PggeRI*GRI-cost_gas_RI)-lambda*(CIg-LCFS)*Sg;
chng_PS2 = chng_PS_corn_LCFS + chng_PS_agres_LCFS + chng_PS_orch_vin_LCFS + chng_PS_forest_LCFS + chng_PS_msw_LCFS + chng_PS_msw_food_LCFS + chng_PS_hec_LCFS + chng_PS_gas_LCFS;

% Save data for combined graphs
Supply_Gas_LCFS = Sg/10^9;         % Convert to billion gge
Supply_EtOH_LCFS = Se/10^9;        % Convert to billion gge
Supply_GGE_LCFS = Sgge/10^9;          % Convert to billion gge
Price_GGE_LCFS = Pgge;
Shadow_Price_LCFS = lambda/(CFg/(10^6));        % Units of $/MT
Supply_EtOH_corn_LCFS = Se_corn/10^9;          % Convert to billion gge
Supply_EtOH_agres_LCFS = Se_agres/10^9;            % Convert to billion gge
Supply_EtOH_orch_vin_LCFS = Se_orch_vin/10^9;          % Convert to billion gge
Supply_EtOH_forest_LCFS = Se_forest/10^9;             % Convert to billion gge
Supply_EtOH_hec_LCFS = Se_hec/10^9;            % Convert to billion gge
Supply_EtOH_msw_LCFS = Se_msw/10^9;            % Convert to billion gge
Supply_EtOH_msw_food_LCFS = Se_msw_food/10^9;            % Convert to billion gge
Baseline_Carbon_LCFS = C_RI/10^6; % Convert to million metric tons
Carbon_LCFS= C/10^6;                  % Convert to million metric tons
Change_CS_LCFS = chng_CS/10^9;           % In billions of $
Change_PS_LCFS = chng_PS/10^9;           % In billions of $
Abatement_Cost_LCFS = chng_SS/10^9;           % In billions of $
Avg_Abatement_cost_LCFS = Avg_abatement_cost;  % In $/MTCO2e

% Save prive vector & CSPSVec
PPriceVec_LCFS = PPriceVec;
CSPSVec_LCFS= [chng_CS;chng_PS_corn_LCFS;chng_PS_agres_LCFS;chng_PS_orch_vin_LCFS;chng_PS_forest_LCFS;chng_PS_msw_LCFS;chng_PS_msw_food_LCFS;chng_PS_hec_LCFS;chng_PS_gas_LCFS;99999999];



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% VEETC and Cellulosic Subsidy Calculations
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp 'VEETC and Cellulosic Subsidy Calculations'

Pgge = PggeRI;  
Sg = GRI;               % Initialize gasoline supply
Se_corn = 0;            
Se_agres = 0;
Se_orch_vin = 0;        % Initialize ethanol supplies
Se_forest = 0;
Se_hec = 0;
Se_msw = 0;
Se_msw_food = 0;
Sgge = GGERI;           % Initialize supply of BTU's

% Calculate cost in BAU
Cost_no_subsidy = Cost_RI;   % Set cost to baseline cost in other scenarios

% Adjust ethanol supply curves for blender's credit
VEETC = 0.45*BTUPGg/BTUPGe; % Ethanol (corn) blender's credit, source = AFDC
Cellulosic_subsidy = 1.01*BTUPGg/BTUPGe; % Ethanol (cellulosic) blender's credit, source = AFDC

% Set fuel price to BAU
Pgge = PggeRI;
PPriceVec=ones(7,1)*(Pgge+Cellulosic_subsidy);
PPriceVec(1)=Pgge+VEETC;
CalculateSupply;

%
  Dgge = GGERI*((Pgge/PggeRI)^Egged);      % Calculate demand for fuel at market price
  Se = Se_corn + Se_agres + Se_orch_vin + Se_forest + Se_hec + Se_msw + Se_msw_food;       % Calculate supply of ethanol
  Sg = Dgge - Se;      % Calculate supply of gasoline given perfectly elastic supply
    if Sg < 0,  disp     'Error supply of gas is negative (CT)'; lambda; LCFS; err5,    end;
        
% Calculate costs
% Carbon in Metric Tons CO2e
C_RI = (GRI*CFg/10^6)+(ERI_corn*CFe_corn/10^6)+(ERI_agres*CFe_agres/10^6)+(ERI_orch_vin*CFe_orch_vin/10^6)+(ERI_forest*CFe_forest/10^6)+(ERI_hec*CFe_hec/10^6)+(ERI_msw*CFe_msw/10^6)+(ERI_msw_food*CFe_msw_food/10^6); % CO2 eq. with LCFS - Result here is MMTCO@e 
C = (Sg*CFg/10^6)+(Se_corn*CFe_corn/10^6)+(Se_agres*CFe_agres/10^6)+(Se_orch_vin*CFe_orch_vin/10^6)+(Se_forest*CFe_forest/10^6)+(Se_hec*CFe_hec/10^6)+(Se_msw*CFe_msw/10^6)+(Se_msw_food*CFe_msw_food/10^6); % CO2 eq. with LCFS - Result here is MMTCO@e

% Calculate change in private surplus and average abatement cost
% Total Cost of Production

    % Gasoline
    cost_gas = PggeRI * Sg;
    
    % Calculate total production costs under policy
    Cost = cost_corn+cost_agres+cost_orch_vin+cost_forest+cost_hec+cost_msw+cost_msw_food+cost_gas;

    % Calculate supply of fuel/GGE's
    Sgge = Sg + Se;

% Area under demand curves
% Set lower bound for integral (so that integral is finite
Dem_area = PggeRI*(1/GGERI)^(1/Egged)*(1/(1/Egged+1))*Dgge^(1/Egged+1) - PggeRI*(1/GGERI)^(1/Egged)*(1/(1/Egged+1))*int_limit^(1/Egged+1);

% Calculations
SS_RI = Dem_area_RI - Cost_RI;
SS = Dem_area - Cost;
chng_SS = SS-SS_RI;
chng_carbon = C-C_RI;
Avg_abatement_cost = chng_SS/chng_carbon;

% Consumer Surplus
CS_RI = PggeRI*(1/GGERI)^(1/Egged)*(1/(1/Egged+1))*GGERI^(1/Egged+1) - PggeRI*(1/GGERI)^(1/Egged)*(1/(1/Egged+1))*int_limit^(1/Egged+1) - (PggeRI*(GGERI-int_limit));
CS =  PggeRI*(1/GGERI)^(1/Egged)*(1/(1/Egged+1))*Dgge^(1/Egged+1) - PggeRI*(1/GGERI)^(1/Egged)*(1/(1/Egged+1))*int_limit^(1/Egged+1) - (Pgge*(Dgge-int_limit));
chng_CS = CS - CS_RI;


% Calculate subsidy payments
Subsidy_Payments = Se_corn*VEETC + Cellulosic_subsidy*(Se_agres + Se_orch_vin + Se_forest + Se_hec + Se_msw + Se_msw_food); 

% Producer Surplus
PS_RI = PggeRI*GGERI - Cost_RI;
PS = Pgge*Sgge - Cost;
% Add in subsidies
PS = PS + Subsidy_Payments;
chng_PS = PS - PS_RI; 

% Producer Surplus by fuel type
chng_PS_corn_SUBS = (Pgge*Se_corn-cost_corn)-(PggeRI*ERI_corn-cost_corn_RI)+VEETC*Se_corn;
chng_PS_agres_SUBS = (Pgge*Se_agres-cost_agres)-(PggeRI*ERI_agres-cost_agres_RI)+Cellulosic_subsidy*Se_agres;
chng_PS_orch_vin_SUBS = (Pgge*Se_orch_vin-cost_orch_vin)-(PggeRI*ERI_orch_vin-cost_orch_vin_RI)+Cellulosic_subsidy*Se_orch_vin;
chng_PS_forest_SUBS = (Pgge*Se_forest-cost_forest)-(PggeRI*ERI_forest-cost_forest_RI)+Cellulosic_subsidy*Se_forest;
chng_PS_msw_SUBS = (Pgge*Se_msw-cost_msw)-(PggeRI*ERI_msw-cost_msw_RI)+Cellulosic_subsidy*Se_msw;
chng_PS_msw_food_SUBS = (Pgge*Se_msw_food-cost_msw_food)-(PggeRI*ERI_msw_food-cost_msw_food_RI)+Cellulosic_subsidy*Se_msw_food;
chng_PS_hec_SUBS = (Pgge*Se_hec-cost_hec)-(PggeRI*ERI_hec-cost_hec_RI)+Cellulosic_subsidy*Se_hec;
chng_PS_gas_SUBS = (Pgge*Sg-cost_gas)-(PggeRI*GRI-cost_gas_RI);
chng_PS2 = chng_PS_corn_SUBS + chng_PS_agres_SUBS + chng_PS_orch_vin_SUBS + chng_PS_forest_SUBS + chng_PS_msw_SUBS + chng_PS_msw_food_SUBS + chng_PS_hec_SUBS + chng_PS_gas_SUBS;

% Save data for combined graphs
Supply_Gas_SUB = Sg/10^9;         % Convert to billion gge
Supply_EtOH_SUB = Se/10^9;        % Convert to billion gge
Supply_GGE_SUB = Sgge/10^9;          % Convert to billion gge
Price_GGE_SUB = Pgge;
Supply_EtOH_corn_SUB = Se_corn/10^9;          % Convert to billion gge
Supply_EtOH_agres_SUB = Se_agres/10^9;            % Convert to billion gge
Supply_EtOH_orch_vin_SUB = Se_orch_vin/10^9;          % Convert to billion gge
Supply_EtOH_forest_SUB = Se_forest/10^9;             % Convert to billion gge
Supply_EtOH_hec_SUB = Se_hec/10^9;            % Convert to billion gge
Supply_EtOH_msw_SUB = Se_msw/10^9;            % Convert to billion gge
Supply_EtOH_msw_food_SUB = Se_msw_food/10^9;            % Convert to billion gge
Baseline_Carbon_SUB = C_RI/10^6; % Convert to million metric tons
Carbon_SUB= C/10^6;                  % Convert to million metric tons
Change_CS_SUB = chng_CS/10^9;           % In billions of $
Change_PS_SUB = chng_PS/10^9;           % In billions of $
Abatement_Cost_SUB = chng_SS/10^9;           % In billions of $
Avg_Abatement_cost_SUB = Avg_abatement_cost;  % In $/MTCO2e
Ethanol_Subsidies = Subsidy_Payments/10^9; % In billions of $

% Save prive vector & CSPSVec
PPriceVec_SUBS = PPriceVec;
CSPSVec_SUBS= [chng_CS;chng_PS_corn_SUBS;chng_PS_agres_SUBS;chng_PS_orch_vin_SUBS;chng_PS_forest_SUBS;chng_PS_msw_SUBS;chng_PS_msw_food_SUBS;chng_PS_hec_SUBS;chng_PS_gas_SUBS;Subsidy_Payments];



% Output Simulation Results    csvwrite(filename,M,row,col)
if RFStoggle==1
    output_rfs = [Supply_Gas_RFS;Supply_EtOH_RFS;Supply_GGE_RFS;Price_GGE_RFS;Supply_EtOH_corn_RFS;Supply_EtOH_agres_RFS;Supply_EtOH_orch_vin_RFS;Supply_EtOH_forest_RFS;Supply_EtOH_hec_RFS;Supply_EtOH_msw_RFS;Supply_EtOH_msw_food_RFS;Baseline_Carbon_RFS;Carbon_RFS;Change_CS_RFS;Change_PS_RFS;Abatement_Cost_RFS;Avg_Abatement_cost_RFS;lambda_corn;lambda_adv;lambda_cell];
else
    output_rfs=zeros(20,1);
end
output_ct = [Supply_Gas_CT;Supply_EtOH_CT;Supply_GGE_CT;Price_GGE_CT;Supply_EtOH_corn_CT;Supply_EtOH_agres_CT;Supply_EtOH_orch_vin_CT;Supply_EtOH_forest_CT;Supply_EtOH_hec_CT;Supply_EtOH_msw_CT;Supply_EtOH_msw_food_CT;Baseline_Carbon_CT;Carbon_CT;Change_CS_CT;Change_PS_CT;Abatement_Cost_CT;Avg_Abatement_cost_CT;Shadow_Price_CT;CT_Revenue_CT];
output_lcfs = [Supply_Gas_LCFS;Supply_EtOH_LCFS;Supply_GGE_LCFS;Price_GGE_LCFS;Supply_EtOH_corn_LCFS;Supply_EtOH_agres_LCFS;Supply_EtOH_orch_vin_LCFS;Supply_EtOH_forest_LCFS;Supply_EtOH_hec_LCFS;Supply_EtOH_msw_LCFS;Supply_EtOH_msw_food_LCFS;Baseline_Carbon_LCFS;Carbon_LCFS;Change_CS_LCFS;Change_PS_LCFS;Abatement_Cost_LCFS;Avg_Abatement_cost_LCFS;Shadow_Price_LCFS];
output_bau = [GRI/10^9;ERI/10^9;GGERI/10^9;PggeRI;ERI_corn/10^9;ERI_agres/10^9;ERI_orch_vin/10^9;ERI_forest/10^9;ERI_hec/10^9;ERI_msw/10^9;ERI_msw_food/10^9;Baseline_Carbon/10^6;];
output_subsidy = [Supply_Gas_SUB;Supply_EtOH_SUB;Supply_GGE_SUB;Price_GGE_SUB;Supply_EtOH_corn_SUB;Supply_EtOH_agres_SUB;Supply_EtOH_orch_vin_SUB;Supply_EtOH_forest_SUB;Supply_EtOH_hec_SUB;Supply_EtOH_msw_SUB;Supply_EtOH_msw_food_SUB;Baseline_Carbon_SUB;Carbon_SUB;Change_CS_SUB;Change_PS_SUB;Abatement_Cost_SUB;Avg_Abatement_cost_SUB;Ethanol_Subsidies];

% Output PS by Fuel Results
if RFStoggle==1
    ps_RFS = [chng_PS_corn_RFS;chng_PS_agres_RFS;chng_PS_orch_vin_RFS;chng_PS_forest_RFS;chng_PS_msw_RFS;chng_PS_msw_food_RFS;chng_PS_hec_RFS;chng_PS_gas_RFS];
else
    ps_RFS=zeros(8,1);
end
ps_CAT = [chng_PS_corn_CT;chng_PS_agres_CT;chng_PS_orch_vin_CT;chng_PS_forest_CT;chng_PS_msw_CT;chng_PS_msw_food_CT;chng_PS_hec_CT;chng_PS_gas_CT];
ps_LCFS = [chng_PS_corn_LCFS;chng_PS_agres_LCFS;chng_PS_orch_vin_LCFS;chng_PS_forest_LCFS;chng_PS_msw_LCFS;chng_PS_msw_food_LCFS;chng_PS_hec_LCFS;chng_PS_gas_LCFS];
ps_SUBS = [chng_PS_corn_SUBS;chng_PS_agres_SUBS;chng_PS_orch_vin_SUBS;chng_PS_forest_SUBS;chng_PS_msw_SUBS;chng_PS_msw_food_SUBS;chng_PS_hec_SUBS;chng_PS_gas_SUBS];


if user==1    

if Innov_toggle == 0

% Create output vectors
output_all = [output_bau;output_rfs;output_lcfs;output_ct;output_subsidy;ps_RFS;ps_LCFS;ps_CAT;ps_SUBS;PS_RI;CS_RI];
PriceVectors = [PPriceVec_BAU PPriceVec_RFS PPriceVec_CT PPriceVec_LCFS PPriceVec_SUBS];
CSPSVectors = [CSPSVec_BAU CSPSVec_RFS CSPSVec_CT CSPSVec_LCFS CSPSVec_SUBS]; 

if scenario == 0
csvwrite('/Users/Jon/Dropbox/LandUse/2 - Policy Simulation/Simulation Results/RFS2SimulationALL_new.csv',output_all,0,0)
csvwrite('/Users/Jon/Dropbox/LandUse/2 - Policy Simulation/Simulation Results/PriceVectors_new.csv',PriceVectors);
csvwrite('/Users/Jon/Dropbox/LandUse/2 - Policy Simulation/Simulation Results/CSPSVectors_new.csv',CSPSVectors);
end
if scenario == 1
csvwrite('/Users/Jon/Dropbox/LandUse/2 - Policy Simulation/Simulation Results/RFS2SimulationALL_lowP.csv',output_all,0,0)
csvwrite('/Users/Jon/Dropbox/LandUse/2 - Policy Simulation/Simulation Results/PriceVectors_lowP.csv',PriceVectors);
csvwrite('/Users/Jon/Dropbox/LandUse/2 - Policy Simulation/Simulation Results/CSPSVectors_lowP.csv',CSPSVectors);
end
if scenario == 2
csvwrite('/Users/Jon/Dropbox/LandUse/2 - Policy Simulation/Simulation Results/RFS2SimulationALL_highP.csv',output_all,0,0)
csvwrite('/Users/Jon/Dropbox/LandUse/2 - Policy Simulation/Simulation Results/PriceVectors_highP.csv',PriceVectors);
csvwrite('/Users/Jon/Dropbox/LandUse/2 - Policy Simulation/Simulation Results/CSPSVectors_highP.csv',CSPSVectors);
end
if scenario == 3
csvwrite('/Users/Jon/Dropbox/LandUse/2 - Policy Simulation/Simulation Results/RFS2SimulationALL_highiLUC.csv',output_all,0,0)
csvwrite('/Users/Jon/Dropbox/LandUse/2 - Policy Simulation/Simulation Results/PriceVectors_highiLUC.csv',PriceVectors);
csvwrite('/Users/Jon/Dropbox/LandUse/2 - Policy Simulation/Simulation Results/CSPSVectors_highiLUC.csv',CSPSVectors);
end
if scenario == 4
csvwrite('/Users/Jon/Dropbox/LandUse/2 - Policy Simulation/Simulation Results/RFS2SimulationALL_WastePass.csv',output_all,0,0)
csvwrite('/Users/Jon/Dropbox/LandUse/2 - Policy Simulation/Simulation Results/PriceVectors_WastePass.csv',PriceVectors);
csvwrite('/Users/Jon/Dropbox/LandUse/2 - Policy Simulation/Simulation Results/CSPSVectors_WastePass.csv',CSPSVectors);
end
if scenario == 5
csvwrite('/Users/Jon/Dropbox/LandUse/2 - Policy Simulation/Simulation Results/RFS2SimulationALL_OldCorn.csv',output_all,0,0)
csvwrite('/Users/Jon/Dropbox/LandUse/2 - Policy Simulation/Simulation Results/PriceVectors_OldCorn.csv',PriceVectors);
csvwrite('/Users/Jon/Dropbox/LandUse/2 - Policy Simulation/Simulation Results/CSPSVectors_OldCorn.csv',CSPSVectors);
end
if scenario == 6
csvwrite('/Users/Jon/Dropbox/LandUse/2 - Policy Simulation/Simulation Results/RFS2SimulationALL_LessElastic.csv',output_all,0,0)
csvwrite('/Users/Jon/Dropbox/LandUse/2 - Policy Simulation/Simulation Results/PriceVectors_LessElastic.csv',PriceVectors);
csvwrite('/Users/Jon/Dropbox/LandUse/2 - Policy Simulation/Simulation Results/CSPSVectors_LessElastic.csv',CSPSVectors);
end
if scenario == 7
csvwrite('/Users/Jon/Dropbox/LandUse/2 - Policy Simulation/Simulation Results/RFS2SimulationALL_MoreElastic.csv',output_all,0,0)
csvwrite('/Users/Jon/Dropbox/LandUse/2 - Policy Simulation/Simulation Results/PriceVectors_MoreElastic.csv',PriceVectors);
csvwrite('/Users/Jon/Dropbox/LandUse/2 - Policy Simulation/Simulation Results/CSPSVectors_MoreElastic.csv',CSPSVectors);
end
if scenario == 8
csvwrite('/Users/Jon/Dropbox/LandUse/2 - Policy Simulation/Simulation Results/RFS2SimulationALL_CornElast.csv',output_all,0,0)
csvwrite('/Users/Jon/Dropbox/LandUse/2 - Policy Simulation/Simulation Results/PriceVectors_CornElast.csv',PriceVectors);
csvwrite('/Users/Jon/Dropbox/LandUse/2 - Policy Simulation/Simulation Results/CSPSVectors_CornElast.csv',CSPSVectors);
end
end

if Innov_toggle == 1
output_all = [output_bau;output_rfs;output_lcfs;output_ct;output_subsidy;ps_RFS;ps_LCFS;ps_CAT;ps_SUBS;PS_RI;CS_RI;PggeRI*ERI_corn-cost_corn_RI];
csvwrite('/Users/Jon/Dropbox/LandUse/2 - Policy Simulation/Simulation Results/RFS2SimulationALL_NoInnov.csv',output_all,0,0)
end

end


if user==2
    %     if RFStoggle==1
    %         csvwrite('C:\AllStephen\My Dropbox\Jon n Chris\Simulation Results\sph\RFS2SimulationRFS.csv',output_rfs,0,0)
    %     end
    %     csvwrite('C:\AllStephen\My Dropbox\Jon n Chris\2 - Policy Simulation\Simulation Results\sph\RFS2SimulationCT.csv',output_ct,0,0)
    %     csvwrite('C:\AllStephen\My Dropbox\Jon n Chris\2 - Policy Simulation\Simulation Results\sph\RFS2SimulationLCFS.csv',output_lcfs,0,0)
    %     csvwrite('C:\AllStephen\My Dropbox\Jon n Chris\2 - Policy Simulation\Simulation Results\sph\RFS2SimulationBAU.csv',output_bau,0,0)
    %     csvwrite('C:\AllStephen\My Dropbox\Jon n Chris\2 - Policy Simulation\Simulation Results\sph\RFS2SimulationSUBS.csv',output_subsidy,0,0)
    output_all = [output_bau;output_rfs;output_lcfs;output_ct;output_subsidy;ps_RFS;ps_LCFS;ps_CAT;ps_SUBS;PS_RI;CS_RI;PggeRI*ERI_corn-cost_corn_RI];
    csvwrite('C:\AllStephen\My Dropbox\Jon n Chris\2 - Policy Simulation\Simulation Results\sph\RFS2SimulationALL.csv',output_all,0,0)
end

%% this is the end
